
#include <util/types.h>
#include <util/list.h>

#include "../../imp_types.h"
#include "../image.h"

#define min( a, b ) ( a < b ? a : b )

typedef struct _edge {
	Ref2 p[2];
	real x;
	real inv_slope;
	int active;
} Edge;

static int _compare_edges(Edge * e1, Edge * e2) {

	int y1, y2;

	y1 = min(e1->p[0].y, e1->p[1].y);
	y2 = min(e2->p[0].y, e2->p[1].y);

	return y1 - y2;
}

static int _compare_x(int * x1, int * x2) {
	return (*x1) - (*x2);
}

real imp_image_roi_weight(ImpMatrix * Gray, Ref2 * p, int n) {

	real * array = Gray->data;
	int width = Gray->columns;

	Edge ET[n]; // Edge table
	int ne = 0; // Number of edges

	int C[n * 4]; // List of segments to fill
	int nc;

	Ref2 p1, p2;
	int i; // Auxiliary index
	int y; // Scan line

	int ymin, ymax; // Range of the scan line

	// = Number of edges
	ne = n;

	// = Create all edges and get the ymax
	int dx;
	int dy;

	ymax = p[0].y; // ymax initialization to find the max y coordinate of the point set
	for (i = 0; i < n; i++) {

		p1 = p[i];
		p2 = p[ (i + 1) == n ? 0 : (i + 1)];

		// = Find ymax
		if (p1.y > ymax) ymax = p1.y;

		// = The first point of each edge has the minimal y coordinate of the two points
		if (p1.y == min(p1.y, p2.y)) {
			ET[i].p[0] = p1;
			ET[i].p[1] = p2;
		} else {
			ET[i].p[0] = p2;
			ET[i].p[1] = p1;
		}

		// = Calculates the edge inv_slope. As the increment is made in vertical direction, the slope must be inverted. slope = dy/dx.
		dx = p1.x - p2.x;
		dy = p1.y - p2.y;
		ET[i].inv_slope = (real) dx / (real) dy;
		ET[i].active = 0;
	}

	// = sort edges in relation to the min y coordinate of its points
	gqsort(ET, ne, sizeof (ET[0]), _compare_edges);

	// = Get ymin and ymax coordinates of edges set. ===========================
	ymin = ET[0].p[0].y;
	// ymax was calculated above


	// = Find intersections with scanline y from ymin to ymax ==================
	int start = 0;
	real sum = 0;

	for (y = ymin; y <= ymax; y++) {

		// = Find the min Edge that intersects with y ==========================
		while (!(ET[start].p[0].y <= y && y <= ET[start].p[1].y)) {
			start++;
		}

		nc = 0;
		// = For each edge get the intersections with active edges =============
		for (i = start; i < ne; i++) {

			if (y == ET[i].p[0].y) { // scanline y intersects with the first point of edge ET[i]

				if (y == ET[i].p[1].y) { // ET[i] is a horizontal line

					ET[i].active = 0;
					ET[i].x = ET[i].p[1].x;
					C[nc] = ET[i].x;
					nc++;

				} else { // ET[i] is a non-horizontal line

					ET[i].active = 1;
					ET[i].x = ET[i].p[0].x;

				}
			}

			if (y == ET[i].p[1].y) { // y intersects with the bigger y coordinate of edge ET[i]

				ET[i].active = 0;
				ET[i].x = ET[i].p[1].x;
				C[nc] = ET[i].x;
				nc++;
			}

			if (ET[i].p[0].y < y && y < ET[i].p[1].y) { // y intersects with ET[i]

				ET[i].x += ET[i].inv_slope; // get the new x using DDA algorithm
				C[nc] = (ET[i].x - (int) ET[i].x) >= 0 ? (int) ET[i].x + .5 : (int) ET[i].x; // round
				nc++;
			}
		}

		// = sort x coordinates in left to right order =========================
		gqsort(C, nc, sizeof (C[0]), _compare_x);

		if (nc % 2 != 0) {
			// This should never happen.
			return 0;
		}

		// = Calculates the weight =============================================
		int _r, idx;
		int x;

		// = For each founded segment C[i] <--> C[i+1] accumulates the value of (x,y)
		for (i = 0, _r = y * width; i < nc; i += 2, _r += width) {

			if ((i + 2) < nc && C[i] == C[i + 2]) continue; // = Remove pair copies

			for (x = C[i]; x < C[i + 1]; x++) {
				idx = x + _r;
				sum += array[idx];
			}
		}
	}
	return sum;
}

List * imp_image_roi_points(ImpMatrix * Gray, Ref2 * p, int n) {

	List * points = list_new();

	real * array = Gray->data;
	int width = Gray->columns;

	Edge ET[n]; // Edge table
	int ne = 0; // Number of edges

	int C[n * 4]; // List of segments to fill
	int nc;

	Ref2 p1, p2;
	int i; // Auxiliary index
	int y; // Scan line

	int ymin, ymax; // Range of the scan line

	// = Number of edges
	ne = n;

	// = Create all edges and get the ymax
	int dx;
	int dy;

	ymax = p[0].y; // ymax initialization to find the max y coordinate of the point set
	for (i = 0; i < n; i++) {

		p1 = p[i];
		p2 = p[ (i + 1) == n ? 0 : (i + 1)];

		// = Find ymax
		if (p1.y > ymax) ymax = p1.y;

		// = The first point of each edge has the minimal y coordinate of the two points
		if (p1.y == min(p1.y, p2.y)) {
			ET[i].p[0] = p1;
			ET[i].p[1] = p2;
		} else {
			ET[i].p[0] = p2;
			ET[i].p[1] = p1;
		}

		// = Calculates the edge inv_slope. As the increment is made in vertical direction, the slope must be inverted. slope = dy/dx.
		dx = p1.x - p2.x;
		dy = p1.y - p2.y;
		ET[i].inv_slope = (real) dx / (real) dy;
		ET[i].active = 0;
	}

	// = sort edges in relation to the min y coordinate of its points
	gqsort(ET, ne, sizeof (ET[0]), _compare_edges);

	// = Get ymin and ymax coordinates of edges set. ===========================
	ymin = ET[0].p[0].y;
	// ymax was calculated above


	// = Find intersections with scanline y from ymin to ymax ==================
	int start = 0;
	real sum = 0;

	for (y = ymin; y <= ymax; y++) {

		// = Find the min Edge that intersects with y ==========================
		while (!(ET[start].p[0].y <= y && y <= ET[start].p[1].y)) {
			start++;
		}

		nc = 0;
		// = For each edge get the intersections with active edges =============
		for (i = start; i < ne; i++) {

			if (y == ET[i].p[0].y) { // scanline y intersects with the first point of edge ET[i]

				if (y == ET[i].p[1].y) { // ET[i] is a horizontal line

					ET[i].active = 0;
					ET[i].x = ET[i].p[1].x;
					C[nc] = ET[i].x;
					nc++;

				} else { // ET[i] is a non-horizontal line

					ET[i].active = 1;
					ET[i].x = ET[i].p[0].x;

				}
			}

			if (y == ET[i].p[1].y) { // y intersects with the bigger y coordinate of edge ET[i]

				ET[i].active = 0;
				ET[i].x = ET[i].p[1].x;
				C[nc] = ET[i].x;
				nc++;
			}

			if (ET[i].p[0].y < y && y < ET[i].p[1].y) { // y intersects with ET[i]

				ET[i].x += ET[i].inv_slope; // get the new x using DDA algorithm
				C[nc] = (ET[i].x - (int) ET[i].x) >= 0 ? (int) ET[i].x + .5 : (int) ET[i].x; // round
				nc++;
			}
		}

		// = sort x coordinates in left to right order =========================
		gqsort(C, nc, sizeof (C[0]), _compare_x);

		if (nc % 2 != 0) {
			// This should never happen.
			return 0;
		}

		// = Calculates the weight =============================================
		Ref2 * ref;
		int x;

		// = For each founded segment C[i] <--> C[i+1] accumulates the value of (x,y)
		for (i = 0; i < nc; i += 2) {

			if ((i + 2) < nc && C[i] == C[i + 2]) continue; // = Remove pair copies

			for (x = C[i]; x <= C[i + 1]; x++) {

				ref = mem_alloc(Ref2);
				ref->x = x;
				ref->y = y;

				list_append(points, ref);
			}
		}
	}
	return points;
}